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^ ' A systematic calculation of magnetization and specific heat contributions 

^, . due to fluctuations of vortex lattice in strongly type II superconductors to 

C^ , precision of 1% is presented. We complete the calculation of the two loop low 

temperature perturbation theory by including the umklapp processes. Then 
the gaussian variational method is adapted to calculation of thermodynamic 
Ph ' characteristics of the 2D and the 3D vortex solids in high magnetic field, 

j^ . Based on it as a starting point for a perturbation theory we calculate the 

leading correction providing simultaneously an estimate of precision. The 
results are compared to existing nonperturbative approaches. 
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I. INTRODUCTION 

Existence of vortex lattice in type II superconductors in magnetic field was predicted by 
Abrikosov and subsequently observed in various types of such superconductors ranging 
from metals to high Tc cuprates. In the original treatment mean field Ginzburg - Landau 
(GL) theory which neglects thermal fluctuations of the vortex matter was used. Thermal 
fluctuations are expected to play much larger role in high T^ superconductors than in the 
low temperature ones because the Ginzburg parameter Gi characterizing fluctuations is 
much larger [|l| . In addition the presence of strong magnetic field and strong anisotropy in 
superconductors like BSCCO effectively reduces their dimensionality thereby further 
enhancing effects of thermal fluctuations. Under these circumstances fluctuations make the 
lattice softer influencing various physical properties like magnetization and specific heat 
and eventually lead to vortex lattice melting into vortex liquid far below the mean field 
phase transition line 0|^ clearly seen in both magnetization |^ and specific heat 
experiments [Q . To develop a theory of these fluctuations even in the case of lowest 
Landau level (LLL) corresponding to regions of the phase diagram "close" to ifc2 is a 
nontrivial task and several different approaches were developed. 

At high temperature end (and thereby in the "vortex liquid" phase) Thouless and Ruggeri 
0J^ proposed a perturbative expansion around a homogeneous (liquid) state in which all 
the "bubble" diagrams are resummed. It was shown in field theory that summation of 
bubble diagrams is equivalent to the gaussian variational approach |^. In this approach 
one searches for a "gaussian" state having the lowest energy. The series provides accurate 
results at high temperatures, but become inapplicable for LLL dimensionless temperature 
ot ~ (T - T^f{H))/{THfl'^ smaller than 2 in 2D and for ar ^ iT - T^f{H))/{THf/^ 
smaller than 1 in 3D both quite far above the melting line. Generally, attempts to extend 
the theory to lower temperature by Pade extrapolation were not successful P). It is in fact 
doubtful whether the perturbative results based on gaussian approximation assuming 
translational invariant liquid state should be attempted at low ay. 

In ref. 0, it was shown that below ot < —5 different gaussian states which are no longer 
translationally invariant have lower energy. We will present the detailed calculation here 
(optimized perturbation theory was used to study for liquid state, see refs. [p|,p!0|). It is in 
general a very nontrivial problem to find an inhomogeneous solutions of the corresponding 
"gap equation" (see section IV). However using previous experience with low temperature 
perturbation theory 0,|l3, the problem can be significantly simplified and solved using 
rapidly convergent "modes" expansion. A consistent perturbation theory should start from 



these states ||T^. We then generalize the approach of ref. |§] by setting up a perturbation 
theory around the gaussian Abrikosov lattice state. 

Magnetization and specific heat contributions due to vortex lattices are calculated in 
perturbation theory around this state to next to leading order. This allows to estimate the 
precision of the calculation. It is the worst, about 1%, near the melting point at a^ = —10 
and becomes better for lower a^- At low temperature the result is consistent with the first 
principles low temperature perturbation theory advanced recently to the two loop order 
||Tl| , [T2[| . The previous two loop calculation is completed by including the umklapp 
processes. One can make several definitive qualitative conclusions using the improved 
accuracy of the results. The LLL scaled specific heat monotonously rises from its mean 



field value of 1//3a at ar = — cc to a slightly higher value of 1.05/Pa where Pa = 1-16 is 
Abrikosov parameter. This is at variance with theory ref. |l^ which uses completely 
different ideas and has a freedom of arbitrarily choosing certain parameters on a 2% 
precision level. Although we calculate the contribution of the LLL only, corrections due to 



higher Landau levels calculated earlier in ref. ||T5|,|T6[ using less sophisticated method can 
be included. 

The paper is organized as follows. The model is defined in section II. In section III a brief 
summary of existing results and umklapp corrected two loop perturbative results in both 
2D and 3D are given. Gaussian approximation and the mode expansion used is described 
in section IV. The basic idea of expansion around the best gaussian state is explained in 
section V. The leading corrections are calculated. Results are presented and compared 
with perturbation theory and other theories in section VI together with conclusions 



II. MODELS 



To describe fluctuations of order parameter in thin films or layered superconductors one 
can start with the Ginzburg-Landau free energy: 
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where A = {By, 0) describes a constant magnetic field (considered nonfiuctuating) in 
Landau gauge and covariant derivative is defined by D = V — i^Aj^o = ^- For strongly 
type II superconductors like the high Tc cuprates {k ~ 100) and not too far from Hc2 (this 
is the range of interest in this paper, for the detailed discussion of the range of 
applicability see |jl5[) magnetic field is homogeneous to a high degree due to superposition 
from many vortices. For simplicity we assume a{T) = aTc{l — t), t = T/Tc, although this 
temperature dependence can be easily modified to better describe the experimental 
Hc2{T). The thickness of a layer is Lz- 

Throughout most of the paper will use the coherence length ^ = \/h?/ (2maf,aTc) as a unit 
of length and — "^j, "' Tc = -^^ as a unit of magnetic field. After the order parameter field 
is rescaled as i})"^ -^ ^^^ip'^, the dimensionless free energy (the Boltzmann factor) is: 
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The dimensionless coefficient describing the strength of fluctuations is 
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where Gi is the Ginzburg number in 2D . When 



i-f-fe 

126 



<< 1, the lowest Landau level 



(LLL) approximation can be used |T^. The model then simplifies due to the LLL 



constraint, —^ip = ^ip, rescaling x -^ x/Vb,y — > y/Vb and {i/jI"^ — ^ \iIj\'^\/tz, one obtains 
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where the 2D LLL reduced temperature 
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is the only parameter in the theory p!7|,|5|]. 

For 3D materials with asymmetry along the z axis the GL model takes a form: 
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number is now given by 
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T^?/) , where 7 = rric/mab is anisotropy. The Ginzburg 
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Within the LLL approximation, and rescaling 



X ^ X 
becomes: 
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The 3D reduced temperature is: 
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From now on we work with rescaled quantities only and related them to measured 
quantities in section V. 

III. PERTURBATION THEORIES AND EXISTING NONPERTURBATIVE 

RESULTS 



A variety of perturbative as well as nonperturbative methods have been used to study this 
seemingly simple model. There are two phases. Neglecting thermal fluctuations, one 
obtains the lowest energy configuration ip = for a-r > and 
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for ax < 0, where oa = \/% is the lattice spacing in our units and [3 a = 1-16. 



A. High temperature expansion in the hquid phase 



Homogeneous "vortex liquid" phase (which is not separated from the "normal" phase by a 
transition) has been studied using high temperature perturbation theory by Thouless and 
Ruggeri 



Unfortunately this asymptotic series (even pushed to a very high orders ||T8 



are applicable only when aj- > 2. In order to extend the results to lower ax attempts have 
been made to Fade resum the series [^ imposing a constraint that the result matches the 
Abrikosov mean field as a^ — >■ — cxd, . However, if no matching to the limit of is imposed, 
the perturbative results cannot be significantly improved ^. Experiments 0J^, Monte 

^ all point out 



Carlo simulations [|T^] and nonperturbative Bragg chain approximation 

that there is a first order melting transition around a^ = —12. If this is the case it is 

difficult to support such a constraint. 



B. Low temperature perturbation theory in the sohd phase. Umklapp processes 

Recently a low temperature perturbation theory around Abrikosov solution eq. (|n|) was 



developed and shown to be consistent up to the two loop order |JTT]JT^,[T6| . Since we will 
use in the present study the same basis and notations and also will compare to the 
perturbative results we recount here few basic expressions. The order parameter field ip is 
divided into a nonfiuctuating (mean field) part and a small fluctuation 
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The field x can be expanded in a basis of quasi - momentum eigen functions on LLL in 2D: 
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Then we diagonalize the quadratic term of free energy eq. (Q) to obtain the spectrum. 
Instead of the complex field x, two "real" fields O and A will be used: 
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where jk = \lk\ exp[z6'fc] and definition of 7^ (and all other definitions of functions) can be 

(15) 



found in Appendix A. The eigenvalues found by Eilenberger in ref. |^T| are 

eA(k) = -ar 
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Pa Pa 



where /3fc,is defined in Appendix A. In particular, when A; ^ [^ 
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The second excitation mode eo has a finite gap. The free energy to the two loop order was 



calculated in |12|, however the umklapp processes were not included in some two loop 
corrections. These processes correspond to momentum nonconserving (up to integer times 
inverse lattice constant) four leg vertices ( see Appendix A eq. (|A2|)). We therefore 
recalculated these coefficients. The result in 2D is (see Fig.la,b): 



/, 



eff2D 



2/3. 



T , \0'T\ 



19.9 



+ c„. 



(17) 



where c^ = (log [(e^(k) eo(k)/a|]) = -2.92. In 3D, 
carried in ref. |Tl[ to umklapp processes) gives: 



similar calculation (extending the one 
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C. Nonperturbative methods 



Few nonperturbative methods have been attempted. Tesanovic and coworkers ||14 



developed a method based on an approximate separation of the two energy scales. The 
larger contribution (98%) is the condensation energy, while the smaller one (2%) describes 
motion of the vortices. The result for energy in 2D is: 
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Corresponding expressions in 3D were also derived. 

No melting phase transition is seen since it belongs to the 2% which cannot be accounted 

for within this approach. There exist several Monte Carlo (MC) simulations of the system 

||T9| . The expression eq. (|^) agrees quite well with high temperature perturbation theory 

and MC simulations and has been used to fit both magnetization and specific heat 

experiments [Q, but only to mean field level agrees with low temperature perturbation 

theory. Expanding eq. (^) in l/a^, one obtains an opposite sign of the one loop 

contribution, see Fig. la and discussion in section VI. 

Other interesting nonperturbative methods include the 1/A^ expansion PB|P^ and the 



"Bragg chain fiuctuation approximation" |20 . 



IV. GAUSSIAN VARIATIONAL APPROACH 



A. General Anzatz 



Gaussian variational approach originated in quantum mechanics and has been developed in 
various forms and areas of physics p5| , ^ . In quantum mechanics it consists of choosing a 



gaussian wave function which has the lowest energy expectation value. When fermionic 
fields are present the approximation corresponds to BCS or Hartree - Fock variational 
state. In scalar field theory one optimizes the quadratic part of the free energy 

- (0" - v") G-^"^ {(j)^ - v^) + V{G, 0") (20) 

= K + V 

To obtain "shift" v*^ and "width of the gaussian" G, one minimizes the gaussian effective 
free energy |^B|], which is an exact upper bound on the energy (see proof in |^3|). The 
result of the gaussian approximation can be thought of as resummation of all the " cacti" 
or loop diagrams HJ^. Further corrections will be obtained in section V by inserting this 
solution for G and taking more terms in the expansion of Z. 

B. 2D Abrikosov vortex lattice 

In our case of one complex field one should consider the most general quadratic form 

K=[ {r{x)-v\x))G-\x,y){i^{y)-v{y)) (21) 
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Assuming hexagonal symmetry (a safe assumption for the present purpose), the shift 
should be proportional to the mean field solution eq. (|T^), v{x) = vip{x), with a constant v 
taken real thanks global U{1) gauge symmetry. On LLL, as in perturbation theory, we will 
use variables 0^ and A^ defined in eq. (Ill) instead of i/j{x) 
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The phase defined after eq. ([14D is quite important for simplification of the problem and 
was introduced for future convenience. The most general quadratic form is 
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with matrix of functions G{k) on Brillouin zone to be determined together with the 
constant v by the variational principle. The corresponding gaussian free energy is 
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where {■■■)f, denotes average over Brillouin zone. The minimization equations are: 
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observation is that after we have inserted the phase exp[— 20^/2] in eq. (p2D using our 
experience with perturbation theory, Gao appears exphcitly only on the right hand side of 
the last equation. It also implicitly appears on the left hand side due to a need to invert 
the matrix G. Obviously GoA{k) = is a solution and in this case the matrix diagonalizes. 
However general solution can be shown to differ from this simple one just by a global 
gauge transformation. Subtracting eq. (^) from eq. ( pT]) and using eq. (0), we observe 
that matrix G~^ has a form: 

Eo{k) EoA{k)\^( E{k) + /\^\-ik\ ^2\lk\ 

EoA{k) EA{k) ) \ A2|7fc| E{k) - A^\-fk\ 

where Ai, A2 are constants. Substituting this into the gaussian energy one finds that it 
depends on Ai, A2 via the combination A = a/A^ + A| only. Therefore without loss of 
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These equations look quite intractable, however they can be simplified. The crucial 
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generality we can set A2 = 0, thereby returning to the Goa = case [|7 



Using this observation the gap equations significantly simplify. The function E{k) and the 
constant A satisfy: 
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Using eq. 



a formula 
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derived in Appendix A and the hexagonal symmetry of the spectrum, one deduces that 
E{k) can be expanded in "modes" 



E{k) = Y,En(3n{k) 



(33) 



The integer n determines the distance of a points on reciprocal lattice from the origin, see 
Fig. 4. and x = 6xp[— a^/2] = exp[— 27r/v^] = 0.0265. One estimates that En ^ X^o-t, 
therefore the coefficients decrease exponentially with n. Note (see Fig. 4) that for some 
integers, for example n = 2, 5, 6, /?„ = 0. Retaining only first s modes will be called "the s 
mode approximation". We miinimized numerically the gaussian energy by varying v, A 
and first few modes of E{k). The sample results for various a^ and number of modes are 
given in Table 1. 

Table 1. 

Mode expansion 2D. 



ax 


1 mode 


2 modes 


3 modes 


-1000 


-446023.8395 


-431171.9948 


-431171.9757 


-300 


-40131.29217 


-38796.0277 


-38796.02297 


-100 


-4450.41636 


-4303.28685 


-4303.28593 


-50 


-1106.51575 


-1070.63806 


-1070.63791 


-20 


-171.678045 


-166.690727 


-166.690827 


-10 


-39.292885 


-38.433571 


-38.433645 


-5 


-7.3153440 


-7.2237197 


-7.2237422 



We see that in the interesting region of not very low temperatures the energy converges 
extremely fast. In practice two modes are quite enough. The results for the gaussian 
energy are plotted on Fig.l and will be compared with other approaches in section VI. 
Furthermore one can show that around ot < —4.6, the gaussian liquid energy is larger than 
the gaussian sohd energy. So naturally when a^ < —4.6, one should use the gaussian solid 
to set up a perturbation theory. For a-p > —4.2, there is no solution for the gap equations. 

C. 3D Abrikosov vortex lattice 



In 3D, we expand in bases of plan waves in the third direction times previously used 
quasi - momentum function: 
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where integration over k is understood as integration over Brillouin zone and over kz- Most 
of the derivation and important observations are intact. The modifications are following 
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The corresponding gaussian free energy density (after integration over k^) is: 
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Minimizing the above energy, gap equations similar to that in 2D can be obtained. One 
finds that 

Eo(k) = E(k) + A|7fc|, 
E^(k) = E(k)-A|7fe|. 

i?(k) can be solved by modes expansion two. We minimized numerically the gaussian 
energy by varying f , A and first few modes of E[k). The sample results of free energy 
density for various ax with 3 modes are given in Table 2. 









Table 2. 

Mode expansion 3D. 








ax 


-300 


-100 


-50 


-30 


-20 


-10 


-5.5 


f 


-38757.2294 


-4283.2287 


-1057.6453 


-372.2690 


-159.5392 


-33.9873 


-6.5103 



In practice two modes are also quite enough in 3D. As in the case of 2D, one can show that 
around a^ < —5.5, the gaussian liquid energy is larger than the gaussian solid energy. So 
naturally when ay < —5.5, one should use the gaussian solid to set up a perturbation 
theory in 3D. When around qt > —5, there is no solution for the gap equations. 
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V. CORRECTIONS TO THE GAUSSIAN APPROXIMATION 

In this section, we calculate the lowest order correction to the gaussian approximation 
(that will be called postgaussian correction), which will determine the precision of the 
gaussian approximation. This is necessary in order to fit experiments and compare with 
low temperature perturbation theory and other nonperturbative methods. 
First we review a general idea behind calculating systematic corrections to the gaussian 
approximation |^. The procedure is rather similar to calculating corrections to the 
Hartree-Fock approximations used in fermionic systems. Gaussian variational principle 
provided us with the best (in a certain sense) quadratic part of the free energy K from 
which the "steepest descent" corrections can be calculated. The free energy is divided into 
the quadratic part an a "small" perturbation V . For a general scalar theory defined in eq. 
(^) it takes a form: 



f = K + aV 
1 

2' 



K = i0«G-i«V'' (38) 



Here the auxiliary parameter a was introduced to set up a perturbation theory. It will be 
set to one at the end of calculation. Expanding logarithm of the statistical sum in powers 
of a 

Z= j'v^''exp{-K)exp{-aV)= f V^"^^ (avY exp{-K), (39) 

•^ -^ n=0 ■ 



one retains only first few terms. It was shown in refs. PB[ that generally only two - particle 



irreducible diagrams contribute to the postgaussian correction. The gaussian 
approximation corresponds to retaining only first two terms, n = 0, 1, while the 
postgaussian correction retains in addition the contribution of order a^. 
Feynman rules in our case are shown on Fig. 5. We have two propagators for fields A and O 
and three and four leg vertices. Using these rules the postgaussian correction is presented 
on Fig. 6 as a set of two and three loop diagrams. The corresponding expressions are given 
in Appendix B. The Brillouin zone averages were computed numerically using the three 
modes gaussian solution of the previous section. Now we turn to discussion of the results. 

VI. RESULTS, COMPARISON WITH OTHER APPROACHES AND 

CONCLUSION. 

Results for LLL scaled energy, magnetization and specific heat in 2D are presented on Fig. 
1, 2 and 3 respectively. 

A. Energy 



The gaussian energy provides a rigorous upper bound on free energy |2^. Fig. la shows the 



2D gaussian energy (the dash - dotted line), which in the range of a^ from —30 to —10 is 

11 



just above the mean field (the sohd gray hne). This is because it correctly accounts for the 
(positive) logarithmic one loop correction of eq. ([l^). In contrast the results of the theory 
by Tesanovic et al [|I^ (the dashed gray line) are lower than the mean field. This reflects 
the fact that although the correct large \aT\ limit is built in, the expansion of the 
expression eq. ( p!9[ ) gives negative coefficient of the log \aT\ term. This is inconsistent with 
both the low temperature perturbation theory and the gaussian approximation. The 
difference between this theory and our result is smaller than 2% only when ax < 30 or at 
small aj- below the 2D melting line (which occurs at ax = —13 according to Monte Carlo 
T^ and phenomenological estimates [|TH|,|I|) where the lines become closer again. It never 



gets larger than 10% though. To effectively quantitatively study the model one has to 
subtract the dominant mean field contribution. This is done in the inset of Fig. la. We 
plot the gaussian result (the dash - dotted line), the one loop perturbative result (the solid 
line) and eq. (|19D (the dashed gray line) in an expanded region —100 < a^ < —10. The 



gaussian approximation is a bit higher than the one loop. 

To determine the precision of the gaussian approximation and compare with the 
perturbative two loop result, we further subtracted the one loop contribution on Fig. lb. 
As expected the postgaussian result is lower than the gaussian though higher than the two 
loop. The difference between the gaussian and the postgaussian approximation in the 
region shown is about |A/| = 0.2, which translates into 0.2% at ax = —30, 0.4% at 
qt = —20 and 2% at ax = —12. The fit for the gauss and postgaussian energy in the 
region —30 < ax < —6 are 

al , , , 19.104 60.5271og|aT| 36.511 

fa2D = -TT^f + 21og laxl + 0.119 ^^^ + — + c, 

^Pa (^t (^x ^t 

al , , , 11.68 60.5271og|aT| 38.705 

U,2D = --^ + 21og \ax\ + 0.068 - -^^^ + — ^ + c.. 

z,pA ax arp arp 

In 3D, similarly one found that 

2,JJA Q-X (^T 



B. Magnetization 

The dimensionless LLL magnetization is defined as 

m{ax) = Jl^q^ (41) 

dax 

and the measure magnetization is 

4.M = -^(|.^P> = -^|.^„|^J-yf, (42) 

cruab cruab 2aTc V 47r 

where ^|) is the order parameter of the original model, and ipr is the rescaled one, which is 
equal to %£(^. Thus 
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A-nM = 7^\^m(aT)- (43) 

We plot the scaled magnetization in region —30 < a^ < —6. Again, the mean field 
contribution dominates, so we subtract it in Fig.2. The solid line is the one loop 
approximation, while the gray line is the two loop approximation. At small negative qt 
the postgaussian (the upper gray dash - dotted line) is very close to the two loop result, 
while the gaussian approximation (the dash - dotted line) is a bit lower. All of these lines 
are above mean field. On the other hand, the result of ref. |jT^ (the gray dashed line) is 



below the mean field. Magnetization jump at the melting point is smaller than our 
precision of 2% at ax = —12. Our result for the gaussian magnetization and the 
postgaussian correction in this range can be conveniently fitted with 



ot _ 2 _ 19.10 133.55 _ 121.05 log[-aT] 

~n 2 3 

Pa ^t ^t ^t 

7.525 59.15 43.64 logf-ar] 



mg2D - 75- - -^ ^ —3 3 

h^A ^T ^T T T 



^^g2D — 2 3 ^ 

Cinn iJjrp Cirp 

respectively. 

Similar discussion for the case of 3D can be deduced from eqs. (|l|) , eq. (^OD and 

e*h h' f huj V'"^ 

where the gaussian scaled magnetization can be obtained by differentiation of eq. (|40|) . We 
didn't attempt to calculate the postgaussian correction in 3D. 



C. Specific heat 

The scaled LLL specific heat is defined as 



d'^feffiar) 
dan 



,(„,) = -ililffi^ (45) 



and the original specific heat is related to the scaled specific heat c in 2D via 



c ' 



Att^^T 



.7Th?abT^ -3t - 1 + 6 , , irri^aTc (-t-l + bf , , 
^' 2m*b'T 2 ^ ' m*b'T 2 ^ ' 



We plot the scaled specific heat divided by the mean field value c^/ = '^/Pa in the range 
—30 < a-r < —6 on Fig. 3. The solid line is one loop approximation, while the gray line is 
the two loop approximation. At large la-rl the postgaussian (gray dashed - dotted line) is 
very closed to the one loop result. Finally the gaussian approximation (dash - dotted line) 
is a bit lower. All these lines are slightly above mean field. On the contrary, the result of 
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ref. [^ (dashed gray line) is below the mean field. Our gaussian result and its correction 
in this range can be conveniently fitted with: 

Cg ^ [2 38.2 521.7 363.2 ln[-ar] 
1 + P A 1 -r + -;3 ZT- + 74 



Acg ^ f 15.05 221.1 130.9 ln[-aT] 

~ PA n h 



1 3 4 4 

C^YYl f \ Ojrp Q/rp Cl/rp 

Qualitatively the gaussian specific heat is consistent with experiments [Q which show that 
the specific heat first raise before dropping sharply beyond the melting point. 

D. Conclusions 

In this paper, we applied the gaussian variational principle to the problem of thermal 
fiuctuations in vortex lattice state. Then the correction to it was calculated perturbatively. 
This generalizes corresponding treatment of fiuctuations in the homogeneous phase (vortex 
liquid) by Thouless and coworkers 0. Also umklapp processes were included in the low 
temperature two loop perturbation theory expression. The results of gaussian perturbative 
and some nonperturbative approaches were compared. We hope that increased sensitivity 
of both magnetization and specific heat experiments will test the precision of the theory. 
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APPENDIX A: 

In this Appendix the basic definitions are collected. Brillouin zone averages of products of 
four quasi - momentum functions are defined by: 

o ^ I 1 2 * -^ 

Pk = < m f-j^f-^ > 

jk = < (<^*)V_-^<^-^ > (Al) 

lk,i = < iplip*_kip _--^ip-^ > . 

We also need a more general product {^*ki'^k2'^*k3Vk4) in order to calculate postgaussian 
corrections. This is just a perturbative four - leg vertex: 



{^l^^k2^t,^ki) =exp 



— - [n^ - rii) + I — nik^y 
2 ^ ' ' oa 



<5'?[ki - k2 + ks - k4]A [ki - k2, k2 - k4] 



A [li, y = Y^ exp[- ' ' ^^' + I (/i, + /2,) Qy - i ihy + hy) Qx] (A2) 

Q 
X exp [i (Zi^ + hx) hy] , 
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where are reciprical lattice vectors: 

Q = midi + m2d2. (A3) 

Here ki — k2 + ks — k4 = nidi + n2d2 is assumed and the basis of reciprocal lattice is 

OA oaVs'^ \Mho-n V. — ko + kq 



d' = 5 !■-* 'd== °.^ 



^. It is dual to the lattice 

ei = {a A, 0), ds = (^, ^^"j . When ki - k2 + kg - k4 7^ nidi + n2d2 the quantity 
vanishes. The delta function differs from the Kroneker: 



Q 



From the above formula, one gets the following expansion of /Sk'- 



(A4) 



X^exp 



-T^ 



Pn{k) 



where X =nidi + n2d2. 

To simplify the minimization equations we used the following general identity. Any sixfold 
{Dq) symmetric function F{k) (namely a function satisfying F{k) = F{k ), where k, k is 
related by a ^ rotation) obeys: 



This can be seen by expanding F in Fourier modes and symmetrizing. 



(A5) 



APPENDIX B: 

In this Appendix we specify Feynman rules and collect expressions for diagrams. The solid 
line Fig. 5a represents O and the dashed line Fig.5b represents A. Fig. 5c is a vertex with 

three O. In the coordinate space, it is 2v (pO (O"*") + c.c. . And Fig. 5c is 
-2iv^+0'^A+ - AivipOO+A+ + c.c. Fig.5g is | \0 (a;)|^ Fig.5h is 200+ {A0+ - 0A+). 
Fig.5i is AOO+AA+ - \o^ {A+f + c.c . 

Other vertices, for example, formulas for diagrams Fig.5e, f, j and k, can be obtained by 
substituting O -^ iA, A -^ —iO from formulas for diagrams Fig.5d, c, h and g respectively. 
The propagator in coordinate space can be written as 



(Bl) 



(0+(x)0+(y)> =An / Eoik)^lix)^%iy)=A7rP+ix,y), 

(O(x)O(i/)) =47r f Eoik)M^)V-kiy)='^7vPoix,y), 
Jk 

(0(x)0+(i/)> =An f Eo{k)^k{x)^l{y) = AnPo{x,y) 
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Functions -P4 (x, y), P^ {x, y), Pa{x, y) can be defined similarly. 

One finds three loops contribution to free energy from two 0000 vertex contraction, see 

Fig-6a, -Ti(i)5/x(/oooo)y: 

/oooo = 4(|Po|'+|P^r + 4|PoP+|'\ . (B2) 



Coordinates are not written explicitly since all of them are the same Po{x,y) etc. 
The contribution from the diagrams Fig. 6b is — ,;^ .5 / {foooa)^ ^^"^^ 

foooa = |Po|' (-16P^P^ + 8P0PI) + (B3) 

1^0 r {-^PoPa + I6P0PI) + c.c. 
The diagrams Fig.6c is -j^^ j^ {fooaa)y and 

fooaa = 16 {\Pof + |P<^r) (|Pa|' + \P1\'') + 

^{[P<:>T[PaT + P'o[P1Y + c.c.) 

-32 {PoP*oP1Pa + c.c.) (B4) 

The diagrams Fig.6f is --^^^ J^ {fooo)y and 



fooo = \Po\ (i6P^<^ (x) ^ (y) + 8P^v (x) V* {y) + c.c.) 

+ |P+|' (SP+y. (x) y. (y) + lQP*oV (x) y.* (y) + C.C.) (B5) 



The diagrams Fig.6h is -^^^ J^ (/ooa)yand 



16{ 

2 



O I ' ^ 



/00a = -8 {Po)'Pl^* (x)^* (y) - 16 (|Po|' + |P^ 

(P+(^ (x) ^ (y) - P^^* (x) v^ (y)) + 8P2P>* (x) ^ (y) 

-32PoPo [Plv^* (x) y. (y) - P>* (x) V* (y)] (B6) 

+C.C 

Other contributions, Fig.6e,d,i,g can be obtained by substituting 

Po ^^ Pa Pi ^^ -PqPa ^-^ -Po in eq. (^, eq. (0), eq. (i|) and eq. (0). 
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Figure captions 



Fig. la 

Scaled free energy of vortex solid. From top to bottom, gaussian approximation (dash 
dotted line), mean field (solid line), theory ref. []T1| (dashed line). Inset: corrections to 
mean field calculated using (from top to bottom) gaussian (dash - dotted line), one loop 
perturbation theory (dotted line) and theory ref. [|T^ (dashed line). 

Fig. lb 

More refined comparison of different approximations to free energy. Mean field as well as 
the one loop perturbative contributions are subtracted. 

Fig. 2 

Thermal fluctuations correction to magnetization of vortex solid. From top to bottom, one 
and two loop perturbation theory (solid lines "pi" and "p2" respectively), gaussian and 
postgaussian approximationa (dash dotted lines "g" and "pg" respectively), theory ref. |T^ 
(dashed line "t"). 

Fig. 3 



Scaled specific heat eq.(P5|) normalized by the mean field. One and two loop perturbation 
theory (solid lines "pi" and "p2" respectively), gaussian and postgaussian approximationa 
(dash dotted lines "g" and "pg" respectively), theory ref. ||14| (dashed line "t"). 

Fig. 4 

Reciprical hexagonal lattice points X belonging to three lowest order "stars" in the mode 
expansion of Pk- 

Fig. 5 

Feynman rules of the low temperature perturbation theory. The solid line (a) denoted the 
O mode propagator, while the dashed line (b) denotes the A mode propagator. Various 
three leg and four leg vertices are presented on (c-f ) and (g-k) respectively. 

Fig. 6 

Contributions to the postgaussian correction to free energy. 
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